GC MAGIC visualization
Marc Elosua-Bayes
abril 22, 2022
Introduction
In this Rmarkdown document we are going to load the MAGIC-denoised data to better visualize genes and ease with the annotation when using specific marker genes. MAGIC was developed by Smita Krishnaswamy’s lab to try to fill in the drop out reads in the spots. MAGIC is a Markov Affinity-based Graph Imputation of Cells used for denoising high-dimensional data most commonly applied to single-cell RNA sequencing data. MAGIC learns the manifold data, using the resultant graph to smooth the features and restore the structure of the data based on their k-nearest neighbors.
Libraries
library(Seurat)
library(dplyr)
library(ggplot2)
library(SPATA2)
library(UCell)
library(stringr)
library(readr)Parameters
set.seed(123)
source(here::here("misc/paths.R"))
source(here::here("utils/bin.R"))
"{gct}/{plt_dir}" %>%
glue::glue() %>%
here::here() %>%
dir.create(path = .,
showWarnings = FALSE,
recursive = TRUE)
"{gct}/{robj_dir}" %>%
glue::glue() %>%
here::here() %>%
dir.create(path = .,
showWarnings = FALSE,
recursive = TRUE)Load data
The data used in this Rmarkdown document comes from 03-clustering_integration.Rmd where the data was integrated.
merged_se <- "{anot}/{robj_dir}/integrated_spatial_annot.rds" %>%
glue::glue() %>%
here::here() %>%
readRDS(file = .)Load MAGIC data from the script MAGIC_denoising-GC.Rmd
magic_df <- "{gct}/{robj_dir}/MAGIC-mtrx.rds" %>%
glue::glue() %>%
here::here() %>%
readRDS(file = .)
# create a new assay to store MAGIC information
magic_assay <- CreateAssayObject(counts = as.matrix(magic_df))
# Subset merged_se to those barcodes used
merged_se <- merged_se[, colnames(magic_df)]
# add this assay to the previously created Seurat object
merged_se[["MAGIC_Spatial"]] <- magic_assayAnalysis
Load gene list
"{gct}/GC_dict.R" %>%
glue::glue() %>%
here::here() %>%
source(file = .)
gene_vec <- rownames(magic_df)
gc_vec <- unique(c(gc_dict[["HARD"]], gc_dict[["MID"]]))Marker gene location
Look at the location where the genes of interest are expressed
Seurat::DefaultAssay(merged_se) <- "MAGIC_Spatial"
# Iterate over each image
lapply(id_sp_df$gem_id, function(i) {
print(i)
gene_plt <- Seurat::SpatialFeaturePlot(
object = merged_se,
features = gc_vec,
alpha = c(0, 1),
ncol = 6,
images = i)
# Save plot
"{gct}/{plt_dir}/magic_GC_markers_{i}.pdf" %>%
glue::glue() %>%
here::here() %>%
cowplot::save_plot(
filename = .,
plot = gene_plt,
base_height = 35,
base_width = 35)
})## [1] "tarwe1_xott6q"## [1] "c28w2r_7jne4i"## [1] "esvq52_nluss5"## [1] "p7hv1g_tjgmyj"## [1] "gcyl7c_cec61b"## [1] "zrt7gl_lhyyar"## [1] "qvwc8t_2vsr67"## [1] "exvyh1_66caqq"## [[1]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/magic_GC_markers_tarwe1_xott6q.pdf"
##
## [[2]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/magic_GC_markers_c28w2r_7jne4i.pdf"
##
## [[3]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/magic_GC_markers_esvq52_nluss5.pdf"
##
## [[4]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/magic_GC_markers_p7hv1g_tjgmyj.pdf"
##
## [[5]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/magic_GC_markers_gcyl7c_cec61b.pdf"
##
## [[6]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/magic_GC_markers_zrt7gl_lhyyar.pdf"
##
## [[7]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/magic_GC_markers_qvwc8t_2vsr67.pdf"
##
## [[8]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/magic_GC_markers_exvyh1_66caqq.pdf"Now with the log-norm expression
Seurat::DefaultAssay(merged_se) <- "Spatial"
lapply(id_sp_df$gem_id, function(i) {
# Iterate over each image
gene_plt <- Seurat::SpatialFeaturePlot(
object = merged_se,
features = gc_vec,
alpha = c(0, 1),
ncol = 6,
images = i)
"{gct}/{plt_dir}/lognorm_GC_markers_{i}.pdf" %>%
glue::glue() %>%
here::here() %>%
cowplot::save_plot(
filename = .,
plot = gene_plt,
base_height = 35,
base_width = 35)
})## [[1]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/lognorm_GC_markers_tarwe1_xott6q.pdf"
##
## [[2]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/lognorm_GC_markers_c28w2r_7jne4i.pdf"
##
## [[3]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/lognorm_GC_markers_esvq52_nluss5.pdf"
##
## [[4]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/lognorm_GC_markers_p7hv1g_tjgmyj.pdf"
##
## [[5]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/lognorm_GC_markers_gcyl7c_cec61b.pdf"
##
## [[6]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/lognorm_GC_markers_zrt7gl_lhyyar.pdf"
##
## [[7]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/lognorm_GC_markers_qvwc8t_2vsr67.pdf"
##
## [[8]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/GC-Analysis/2020-09-22/plots_2020-09-22/lognorm_GC_markers_exvyh1_66caqq.pdf"se_sub <- subset(merged_se, subset = gem_id == "esvq52_nluss5")
se_sub## An object of class Seurat
## 29266 features across 3079 samples within 2 assays
## Active assay: Spatial (26846 features, 2000 variable features)
## 1 other assay present: MAGIC_Spatial
## 3 dimensional reductions calculated: pca, umap, harmonyse_sub@images <- se_sub@images[Seurat::Images(se_sub) == "esvq52_nluss5"]Correlation matrix
Since we are working with sample esvq52_nluss5 in this example we will limit the correlation plot to this slide.
(cor_mtrx <- SCrafty::correlation_heatmap(
se = se_sub,
genes = gc_vec,
assay = "MAGIC_Spatial",
slot = "data"))
"{gct}/{plt_dir}/magic_cor-mtrx_markers.pdf" %>%
glue::glue() %>%
here::here() %>%
cowplot::save_plot(
filename = .,
plot = cor_mtrx,
base_height = 15,
base_width = 15)
# Correlation with lognorm expression
cor_log <- SCrafty::correlation_heatmap(
se = se_sub,
genes = gc_vec,
assay = "Spatial",
slot = "data")
"{gct}/{plt_dir}/lognorm_cor-mtrx_markers.pdf" %>%
glue::glue() %>%
here::here() %>%
cowplot::save_plot(
filename = .,
plot = cor_log,
base_height = 9,
base_width = 10)Look at them side by side
cor_mtrx + cor_log
Gene signatures
se_sub <- AddModuleScore_UCell(
obj = se_sub,
features = gc_dict[c("HARD", "MID", "all-targets-intersect")],
name = '_UCell'
)Look at signatures
sign <- colnames(se_sub@meta.data)[stringr::str_detect(
string = colnames(se_sub@meta.data),
pattern = '_UCell')]
Seurat::SpatialPlot(
object = se_sub,
features = sign,
alpha = c(0, 1),
images = "esvq52_nluss5",
ncol = 3)
Follicle focus
follicle_bc <- "{gct}/{robj_dir}/follicle_coordinates.csv" %>%
glue::glue() %>%
here::here() %>%
read_csv()
interzone_bc <- "{gct}/{robj_dir}/LZ_DZ_interzone.csv" %>%
glue::glue() %>%
here::here() %>%
read_csv()
se_follicle <- se_sub[, follicle_bc$barcode]
se_follicle$interzone <- colnames(se_follicle) %in% interzone_bc$barcodeLook at the signature witin the follicle
Seurat::SpatialPlot(
object = se_follicle,
features = sign,
alpha = c(0, 1),
images = "esvq52_nluss5",
ncol = 2,
pt.size.factor = 5) +
Seurat::SpatialPlot(
object = se_follicle,
group.by = "interzone",
images = "esvq52_nluss5",
pt.size.factor = 5)
DE between Interzone & Follicle
Seurat::Idents(object = se_follicle) <- se_follicle@meta.data[, "interzone"]
markers_follicle <- Seurat::FindAllMarkers(
object = se_follicle,
assay = "Spatial",
slot = "data",
verbose = TRUE,
only.pos = TRUE)
# Look at cluster = TRUE
DT::datatable(markers_follicle, filter = "top")Intersection with lists of interest
iz_mgs <- markers_follicle %>% filter(cluster == TRUE) %>% pull(gene)
intersect(iz_mgs, gc_dict[["HARD"]])## character(0)intersect(iz_mgs, gc_dict[["MID"]])## character(0)intersect(iz_mgs, gc_dict[["all-targets-intersect"]])## character(0)Session Info
sessionInfo()## R version 4.0.1 (2020-06-06)
## Platform: x86_64-conda_cos6-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux Server release 6.7 (Santiago)
##
## Matrix products: default
## BLAS/LAPACK: /scratch/groups/singlecell/software/anaconda3/envs/spatial_r/lib/libopenblasp-r0.3.12.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=es_ES.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=es_ES.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=es_ES.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] readr_2.1.2 stringr_1.4.0 UCell_1.2.0 SPATA2_0.1.0 ggplot2_3.3.5 dplyr_1.0.8 SeuratObject_4.0.4 Seurat_4.1.0
##
## loaded via a namespace (and not attached):
## [1] corrplot_0.92 systemfonts_1.0.4 plyr_1.8.6 igraph_1.2.11 lazyeval_0.2.2 splines_4.0.1 crosstalk_1.2.0 listenv_0.8.0 scattermore_0.8 GenomeInfoDb_1.26.7 digest_0.6.29 htmltools_0.5.2 fansi_1.0.2 magrittr_2.0.2 tensor_1.5 cluster_2.1.2 ROCR_1.0-11 limma_3.46.0 tzdb_0.2.0 globals_0.14.0 matrixStats_0.61.0 vroom_1.5.7 spatstat.sparse_2.1-0 colorspace_2.0-3 ggrepel_0.9.1 textshaping_0.3.6 xfun_0.30 crayon_1.5.0 RCurl_1.98-1.6 jsonlite_1.8.0 spatstat.data_2.1-2 survival_3.3-1 zoo_1.8-9 glue_1.6.2 polyclip_1.10-0 gtable_0.3.0 zlibbioc_1.36.0 XVector_0.30.0 leiden_0.3.9 DelayedArray_0.16.3 future.apply_1.8.1 SingleCellExperiment_1.12.0
## [43] BiocGenerics_0.36.1 abind_1.4-5 scales_1.1.1 DBI_1.1.2 spatstat.random_2.1-0 miniUI_0.1.1.1 Rcpp_1.0.8 viridisLite_0.4.0 xtable_1.8-4 reticulate_1.24 spatstat.core_2.4-0 bit_4.0.4 DT_0.21 stats4_4.0.1 htmlwidgets_1.5.4 httr_1.4.2 RColorBrewer_1.1-2 ellipsis_0.3.2 ica_1.0-2 farver_2.1.0 pkgconfig_2.0.3 sass_0.4.0 uwot_0.1.11 deldir_1.0-6 here_1.0.1 utf8_1.2.2 ggcorrplot_0.1.3 labeling_0.4.2 tidyselect_1.1.2 rlang_1.0.2 reshape2_1.4.4 later_1.3.0 munsell_0.5.0 tools_4.0.1 cli_3.2.0 generics_0.1.2 ggridges_0.5.3 evaluate_0.15 fastmap_1.1.0 ragg_1.2.2 yaml_2.3.5 goftest_1.2-3
## [85] bit64_4.0.5 knitr_1.37 fitdistrplus_1.1-6 purrr_0.3.4 RANN_2.6.1 pbapply_1.5-0 future_1.24.0 nlme_3.1-155 mime_0.12 SCrafty_0.1.0 compiler_4.0.1 plotly_4.10.0 png_0.1-7 spatstat.utils_2.3-0 tibble_3.1.6 bslib_0.3.1 stringi_1.7.6 highr_0.9 lattice_0.20-45 Matrix_1.4-0 vctrs_0.3.8 pillar_1.7.0 lifecycle_1.0.1 spatstat.geom_2.3-2 lmtest_0.9-39 jquerylib_0.1.4 RcppAnnoy_0.0.19 data.table_1.14.2 cowplot_1.1.1 bitops_1.0-7 irlba_2.3.5 httpuv_1.6.5 patchwork_1.1.1 GenomicRanges_1.42.0 R6_2.5.1 promises_1.2.0.1 KernSmooth_2.23-20 gridExtra_2.3 IRanges_2.24.1 parallelly_1.30.0 codetools_0.2-18 MASS_7.3-55
## [127] assertthat_0.2.1 SummarizedExperiment_1.20.0 rprojroot_2.0.2 withr_2.5.0 sctransform_0.3.3 S4Vectors_0.28.1 GenomeInfoDbData_1.2.4 mgcv_1.8-39 parallel_4.0.1 hms_1.1.1 grid_4.0.1 rpart_4.1.16 tidyr_1.2.0 rmarkdown_2.12 MatrixGenerics_1.2.1 Rtsne_0.15 Biobase_2.50.0 shiny_1.7.1